

In [2]:

    
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import netCDF4
from datetime import datetime, timedelta

Load the list of files



In [3]:

    
datapath = '/net/atmos/data/erainterim/cdf/2011/10/'
start = datetime(2011,10,5,0)
end = datetime(2011,10,10,12)
def interval(start, end, delta):
    curr = start
    while curr <= end:
        yield curr
        curr += delta
files = [ datapath + date.strftime('Z%Y%m%d_%H') for date in interval(start, end, timedelta(hours=6))]



In [4]:

    
files









    Out[4]:





['/net/atmos/data/erainterim/cdf/2011/10/Z20111005_00',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111005_06',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111005_12',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111005_18',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111006_00',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111006_06',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111006_12',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111006_18',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111007_00',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111007_06',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111007_12',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111007_18',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111008_00',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111008_06',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111008_12',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111008_18',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111009_00',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111009_06',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111009_12',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111009_18',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111010_00',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111010_06',
 '/net/atmos/data/erainterim/cdf/2011/10/Z20111010_12']

Read the files

Use the netCDF4.MFDataset to read a list of file
Use the aklay variable in the pressure_cst file as pressure



In [4]:

    
pfile = netCDF4.Dataset(datapath + 'pressure_cst')
print pfile.variables.keys()
pressure = pfile.variables['aklay'][:]
print pressure
pfile.close()









    



[u'pollon', u'pollat', u'aklev', u'bklev', u'aklay', u'bklay', u'lonmin', u'lonmax', u'latmin', u'latmax', u'dellon', u'dellat', u'starty', u'startm', u'startd', u'starth', u'starts', u'dattyp', u'datver', u'cstver']
[ 900.  850.  800.  700.  600.  500.  400.  300.  250.  200.  100.]



In [5]:

    
ncfile = netCDF4.MFDataset(files)

Load the data

load temperature and specific humidity



In [6]:

    
t = ncfile.variables['T'][:]
q = ncfile.variables['Q'][:]

Compute the equivalent potential temperature using the following information:

Calculate it for the 850 hPa pressure level

Simple equation of the potential temperature $\Theta_{e}$: \begin{equation} \Theta_{e} = T \cdot (1000/P)^{0.286} + 3\cdot w \end{equation}

with T, temperature in K, P the pressure in hPa, w, mixing ratio in g kg$^{-1}$.

The mixing ratio (kg kg$^{-1}$) can be expressed as : \begin{equation} w = \frac{1}{\frac{1}{q_{v}}-1} \end{equation} with $q_{v}$, the specific humidity in kg kg $^{-1}$



In [7]:

    
p = 850.
index = np.where(pressure == p)[0][0]
the = (t[:, index, ...] + 273)*(1000/p)**(0.286) + 3 *1e3 *(1/((1/q[:, index, ...])-1))

Plot $\Theta_{e}$ on a map

Use the domain 50W / 20N to 50E / 80N
Use m.contourf
np.meshgrid can be used to create a grid



In [8]:

    
m = Basemap(resolution='i', llcrnrlon=-50.,llcrnrlat=20.,urcrnrlon=50.,urcrnrlat=80)



In [9]:

    
lon = np.arange(-180,181,1)
lat = np.arange(-90,91,1)
lons, lats = np.meshgrid(lon, lat)



In [10]:

    
m.drawmeridians(np.arange(-180,181,10), labels=[0,0,0,1])
m.drawparallels(np.arange(-90,91,10), labels=[1,0,0,0])
m.drawcoastlines()
m.drawmapboundary()
m.drawcountries()
cf = m.contourf(lons, lats, the[-4, ...], levels = np.arange(285, 335, 5), extend='both', cmap='YlGnBu')
m.colorbar(cf)









    Out[10]:





<matplotlib.colorbar.Colorbar instance at 0x7649368>

Calculate the integrated water vapor flux

$$IVT = \sqrt{\frac{1}{g}\left(\int q_{v}\cdot u\,dp \right)^{2}+ \left(\int q_{v}\cdot v\,dp \right)^{2}}$$

use np.trapz to integrate over pressure level



In [11]:

    
u = ncfile.variables['U'][:]
v = ncfile.variables['V'][:]



In [12]:

    
qu = q * u
qv = q * v



In [13]:

    
quin = 100/(9.81) * np.trapz(qu, pressure, axis=1)
qvin = 100/(9.81) * np.trapz(qv, pressure, axis=1)
ivt = np.sqrt(quin**2 + qvin**2)

Add IVT contours to the $\Theta_{e}$ plot

use m.contour
using the YlGnBu colormap for $\Theta_{e}$ is quite nice



In [14]:

    
m.drawmeridians(np.arange(-180,181,10), labels=[0,0,0,1])
m.drawparallels(np.arange(-90,91,10), labels=[1,0,0,0])
m.drawcoastlines()
m.drawmapboundary()
m.drawcountries()
cf = m.contourf(lons, lats, the[-1, ...], levels = np.arange(285, 335, 5), extend='both', cmap='YlGnBu')
m.colorbar(cf)
cf = m.contour(lons, lats, ivt[-1, ...], colors='r', levels=np.arange(0,700,200))
lb = clabel(cf, cf.levels, fmt='%1.0f')